%% process_field
function process_field(mesh, U)
x = mesh.node(:,1);
y = mesh.node(:,2);
ux = U(1:mesh.Nb);
uy = U(mesh.Nb+1:2*mesh.Nb);
p = U(2*mesh.Nb+1:end);
% interpolation U
xg = linspace(min(x), max(x), 64);
yg = linspace(min(y), max(y), 64);
[X, Y] = meshgrid(xg, yg);
Ux = griddata(x, y, ux, X, Y);
Uy = griddata(x, y, uy, X, Y);
% interpolation P
xp = mesh.node(1:mesh.Np,1);
yp = mesh.node(1:mesh.Np,2);
P = griddata(xp, yp, p, X, Y);
% vector field
figure;
velocityMagnitude = sqrt(Ux.^2 + Uy.^2);
pcolor(X, Y, velocityMagnitude);
shading interp;
colorbar;
hold on;
step = 2;
h = quiver(X(1:step:end, 1:step:end), Y(1:step:end, 1:step:end), Ux(1:step:end, 1:step:end), Uy(1:step:end, 1:step:end), 2);
h.Color = "w";
title("Velocity Field");
hold off;
% streamlines
figure;
streamslice(X, Y, Ux, Uy, 2);
title("Streamlines");
% contours
figure;
contourf(X, Y, P, 20);
colorbar;
title("Pressure Contours");
end